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Fundamental solitons pinned to the interface between two discrete lattices coupled at a single site 
are investigated. Serially and parallel-coupled identical chains (System 1 and System 2), with the 
self-attractive on-site cubic nonlinearity, are considered in one dimension. In these two systems, 
which can be readily implemented as arrays of nonlinear optical waveguides, symmetric, antisym- 
metric and asymmetric solitons are investigated by means of the variational approximation (VA) 
and numerical methods. The VA demonstrates that the antisymmetric solitons exist in the entire 
parameter space, while the symmetric and asymmetric modes can be found below some critical 
value of the coupling parameter. Numerical results confirm these predictions for the symmetric and 
asymmetric fundamental modes. The existence region of numerically found antisymmetric solitons 
is also limited by a certain value of the coupling parameter. The symmetric solitons are destabilized 
via a supercritical symmetry-breaking pitchfork bifurcation, which gives rise to stable asymmetric 
solitons, in both systems. The antisymmetric fundamental solitons, which may be stable or not, do 
, not undergo any bifurcation. In bistability regions stable antisymmetric solitons coexist with either 

£S| , symmetric or asymmetric ones. 
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PACS numbers: 03.75.Lm; 05.45.Yv 



I. INTRODUCTION 



The study of surface modes in multi-layered optical media has begun long ago [l[ . Recently, a great deal of attention 
was attracted to studies of solitons pinned to interfaces between different nonlinear optical media, at least one of which 
carries a lattice structure. Such surface solitons were predicted theoretically in diverse settings [2j and soon after that 
created in experiments Surface solitons were also considered in the general context of junctions between different 
discrete lattices In these studies, it was concluded that the self-trapped surface modes acquire novel properties, 
different from those of the solitons known in uniform lattices. In particular, discrete surface states can only exist 
above a certain threshold value of the total power (the soliton's norm), and a bistability is possible, with different 
■ surface modes coexisting at a common value of the power. 

Generally speaking, these nonlinear modes may be considered as a variety of optical solitons pinned by defects, 
which were also studied theoretically [H, Q and experimentally 0, Q in many systems, continual and discrete (a 
comprehensive review of the topic of discrete and lattice solitons in optics, including the interaction with defects, was 
given in recent articles @; for a general review of discrete solitons, see book [l(j). These studies have demonstrated 
that discrete nonlinear photonic systems may support spatially localized states with different symmetries, which can 
be controlled by the insertion of suitable defects into the lattice @ . Localized modes supported by defects of optical 
lattices (OLs) were also studied in models of Bose-Einstein condensates (BECs) 

Another topic which is relevant to the present work is the possibility of the spontaneous symmetry breaking (SSB) 
in two-mode symmetric settings with a linear coupling between two subsystems. An alternative realization of settings 
which give rise to the SSB is represented by double-well potentials (in that case, the two coupled modes corre- 
spond to states trapped in the two symmetric wells). A specific version of the SSB corresponds to the double- well 
pseudopotential, induced by a symmetric spatial modulation of the local nonlinearity coefficient fl2l |. 

SSB bifurcations, which destabilize symmetric states and give rise to asymmetric ones, were originally predicted in 
terms of the self-trapping in discrete systems [l3| . In the physically important model of dual-core nonlinear fibers, the 
SSB instability was discovered in Ref. [IH, and the respective bifurcations for continuous- wave states were studied 
in detail in Ref. (TBj . for various types of the intra-core nonlinearities. It is also relevant to mention early work [l6j . 
which put forward the SSB concept in the framework of the nonlinear Schrodinger (NLS) equation. Furt her, the SSB 
was studied in detail for solitons in the model of the dual-core fiber with the cubic (Kerr) nonlinearity [l?], EH ■ Similar 
analysis was later performed for gap solitons in the models of dual-core and tri-core [2(| fiber Bragg gratings, 
and for matter-wave solitons in the BEC loaded into a dual-core potential trap combined with a longitudinal OL. 
The latter analysis was performed in the models with both one [2l| and two [22| longitudinal dimensions (ID and 2D, 
respectively) . 

The limit case of a very strong OL corresponds to a discrete lattice [23|. In that case, the SSB of ID and 2D discrete 
solitons in the system of two linearly coupled discrete NLS equations (DNLSEs) was investigated in Ref. [24[ ■ 
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In addition to the studies of the SSB in diverse two-core systems with the cubic nonlinear terms, this effect was also 
analyzed in models describing optical media with quadratic [25j | and cubic-quintic (26j nonlinearities. In the latter 
case, the SSB diagrams feature loops, with asymmetric solitons existing at intermediate values of the total power, 
while only symmetric modes can be found at low and high powers. 

General conclusions about the character of the SSB in solitons can be drawn from the above-mentioned works. 
The self-focusing nonlinearity induces the SSB of symmetric modes, with a trend to make the respective bifurcation 
subcritical, whose characteristic feature is a bistability region, in which stable symmetric and asymmetric solitons 
coexist. The self-focusing nonlinearity does not induce any bifurcation of antisymmetric solitons. In the absence of 
the periodic potential induced by the OL, the entire family of the antisymmetric solitons is unstable. 

However, if the self- attractive nonlinearity acts in the combination with a sufficiently strong OL, the antisymmetric 
solitons may be stable. In the same case, the SSB bifurcation is transformed from subcritical into the supercritical one, 
which does not admit the coexistence of stable symmetric and asymmetric solitons. Nevertheless, a global bistability 
of a different type takes place in this situation, as stable antisymmetric solitons coexist with their symmetric and 
asymmetric counterparts, below and above the bifurcation point, respectively (in the 2D setting, both families of the 
asymmetric and antisymmetric solitons attain a termination point at high values of the total power, due to the onset 
of the collapse driven by the self-attraction). 

On the other hand, in the models combining the self-defocusing nonlinearity and OL potentials, symmetric soli- 
tons do not suffer any SSB, while antisymmetric solitons are destabilized by the anti-symmetry-breaking pitchfork 
bifurcation, which is always of the supercritical type. 

In this work we aim to study localized modes at the interface of two ID uniform discrete lattices (chains) with 
the cubic onsite nonlinearity, which are linearly coupled either in series (System 1, see Fig. Q] below) or in parallel 
(System 2, see Fig. 2]). In fact, the former configuration may be realized as a single discrete lattice with a spring 
defect, in the form of a locally modified inter-site coupling constant (such a system has been actually created in optics, 
as an array of nonlinear waveguides with a modified separation between two of them Q). System 2 differs from that 
introduced in Ref. [24j by the fact that the two uniform identical chains placed in parallel planes are linearly coupled 
in the transverse direction at a singe site, rather than featuring the uniform linear coupling. Analyzing localized 
defect/surface modes in these two systems, we focus on their symmetry properties - in particular, with the intention 
to investigate SSB transitions in them. 

The presentation is organized as follows. The systems, 1 and 2, are introduced in Section 2, where we con- 
sider the existence and stability of various localized modes, applying the variational approximation (VA) and the 
Vakhitov-Kolokolov (VK) stability criterion. Dynamical properties of the discrete solitons with different symmetries 
are investigated by means of numerical methods, including the computation of the linear-stability eigenvalues and 
direct simulations in Section 4. The paper is concluded by Section 5. 



II. FUNDAMENTAL SURFACE SOLITONS IN THE SYSTEM OF TWO COUPLED LATTICES 



A. System 1: the serial coupling 

The system formed by two linked identical semi-infinite chains is displayed in Fig. [TJ with inter-site coupling 
constant C inside the chains, and constant e accounting for the linkage between them. We consider only the case 
when C and s have the same sign, as opposite signs of the inter-site couplings are difficult to realize in optical and 
BEC systems. The on-site nonlinearity is cubic, with the respective coefficient, 7, being constant throughout the 
system. Thus, this model is described by the following DNLSE system, 

i^^ + 4>n+l + 4>n-l + \4>n\ 2 4>n = 0, U ^ 0, 1, 
i^p- + £0x +</>_! + |0 O | 2 00 = 0, 

dz 

i^i+^o + ^ + l^l 2 ^ = 0, (I) 
dz 

where z is the propagation distance (assuming that the chains correspond to two semi-infinite arrays of parallel optical 
waveguides), and an obvious rescaling is used to fix C = 7 = I [that is, e which figures in Eqs. ([!]) is e/C, in terms 

of Fig. m. 

This system may be considered as the usual DNLSE containing a local "spring defect". As mentioned above, such 
a system has been created in optics, in the form of a regular array of parallel waveguiding cores, with a change of the 
distance between two of them [8[ . The interaction of lattice solitons with this array defect was studied experimentally 
in Ref. B. 
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FIG. 1: A schematic of System 1, which consists of two semi-infinite identical chains, linked by the modified linear coupling e 
between sites n — and n — 1. 



Stationary solutions for solitons formed at the interface between the chains are looked for as usual, <p n — u n exp (i/J.z), 
where u n and fi are real lattice field and propagation constant. The corresponding stationary equations following 
from Eqs. ([TJ are 



uu„ + u n+ i + u n -i + = 0, n^O, 1, 
— fiuo + em + u-i + Uq = 0, 
— llui + eu + u 2 + ul = 0. (2) 



1. The variational approximation 

Equations can be derived from the Lagrangian, 

L = L- + L + + 2euqUi, (3) 



n— — oo ^ / \ / 



where L± are the intrinsic Lagrangians of the two semi-infinite chains, and the last term in Eq. ([3]) accounts for 
the coupling between them. To apply the variational approximation (VA), we follow Ref. [2~ij and adopt an ansatz 
consisting of two parts: 

u n = A exp (an), for n < 0, 

u n = B exp (— a(n — 1)), for n > 1. (4) 

This ansatz admits different amplitudes, A =/= £>, of the solution in the linked chains, but postulates a common width 
of the ansatz in both of them, a -1 . Actually, the SSB is represented by the appearance of asymmetric solutions, with 
A 2 ^B 2 . 

Amplitudes A and B are treated below as variational parameters. As concerns inverse width a, it is fixed through 
a solution of the linearization of Eqs. @ at |n| — > oo, 

a = In (,i/2 + vV/4 - l) , (5) 

which implies that the propagation constant takes values fi > 2. Relation ([5]) may also be represented in the following 
forms, that will be used below: 



aS e-« = |-^-l, t i = s + s-\ (6) 

The substitution of ansatz ((4]) into Eq. ([3]) yields the corresponding effective Lagrangian, where Eq. ((5|) is used to 
eliminate fi in favor of s: 

L cti = (L_) cff + (L + ) cS + 2eAB, (7) 

(L)„ = -s- 1 A 2 + , 1 ; A 4 
y ,qS 2 (1 - s 4 ) 

(L + ) cS = -s~ 1 B 2 + 2(1 ^ g 4 ) g4 - (8) 
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This Lagrangian gives rise to the Euler-Lagrange equations for amplitudes A and B: (d/dA) (L-) cS + 2eB = 
0, (d/dB) (L + ) cg + 2eA = 0, or, in the explicit form, 

- s~ 1 A+ - 1 , A 3 +eB = 
1 - s 4 

~ S - 1 B + —^—rB 3 + eA = 0. (9) 

1 - s 4 

These equations allow us to predict the existence of three different species of the localized modes: symmetric ones, 
with A = B, antisymmetric with A — —B, and asymmetric with A 2 ^ B 2 . 

The results of the VA are reported below, along with the respective numerical results. It will be seen that the VA 
predicts characteristics of the symmetric and asymmetric modes very accurately, while the approximation is essentially 
less accurate for the antisymmetric solitons. 



2. Existence regions for the interface solitons 

The solution for the symmetric solitons (SyS) is easily obtained from Eqs. (J5J, 

A 2 = (l-s 4 )(s- 1 -e). (10) 

The dependence of this amplitude on e for fixed s — (5 — v2l) /2 ps 0.21, which corresponds to (J, = 5, is plotted 
in Fig. [21 As follows from Eq. (JTUJ), in System 1 the existence domain of the SyS solutions at a given value of the 
propagation constant, /x, is limited to e <£e = s~~ 1 . The upper existence limit for the SyS is seen in Fig. [21 and also 
in Fig. [3Ja), which shows the existence regions predicted by the VA for all the types of the solitons in the plane of 
(e, //), along with the numerically found counterparts of these regions (the procedure of obtaining numerical results 
is described below). 
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FIG. 2: (Color online) Amplitudes A and B of asymmetric, symmetric, and antisymmetric solitons (black, red, and blue colors, 
respectively), as predicted by the VA (lines) and obtained in the numerical form (AS - circles, SyS - triangles, AnS - squares), 
versus e, for System 1, at fixed values of the propagation constant: [t, = 3 (a) and /x = 5 (b). The dotted green vertical lines 
denote the critical values of e limiting the existence regions of the asymmetric (e c ) and symmetric (e e ) solitons. For fj, = 3, 
antisymmetric solitons were not found in the numerical form. Empty and full symbols pertain to unstable and stable solitons, 
respectively, as concluded from the numerical investigation. 

The existence range of the asymmetric solution (AS), with A 2 ^ B 2 , can be estimated by taking the difference of 
two equations © and canceling a common factor (A — B): 

( s - 1 +e)=A 2 + B 2 +AB. (11) 

In the limit of A — B — > 0, Eq. (flT]) yields the following relation: 

^(l/SXl-s 4 )^- 1 -^. (12) 

Equating expressions (|10[) and (|12[) . one can predict the value of the linkage constant at which the SSB bifurcation 
gives birth to the pair of AS in System 1, 

4 1 ) = (2 S )- 1 . (13) 
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FIG. 3: (Color online) (a) The existence and stability diagrams for the fundamental symmetric (SyS) and asymmetric (AS) 
solitons in System 1, generated by means of the VA (curves) and numerical computations (symbols), in the plane of (e, fi). Black 
squares, which mark the boundary of the SyS existence region, are located very close to the respective line ei 1 (fj) predicted 
by the VA (see Fig. [2]). Red circles, which denote the numerically found boundary of the AS existence region, are also close to 
the VA-predicted curve £^(fx), along which the symmetry-breaking bifurcation takes place, destabilizing the SyS and giving 
birth to two mutually symmetric branches of the AS modes, (b) The existence and stability borders for the antisymmetric 
mode (AnS) in the same parameter plane, found in the numerical form only (the variational approximation does not predict 
existence limits for the AnS). Inside the instability areas in (b), the type of the respective unstable eigenvalues (EVs) of small 
perturbations is indicated. 



For broad solitons, with a — > oo, i.e., s — > 1 [see Eq. (J4j)] , Sc given by Eq. (fl"3j) approaches a limit value, £c — 1/2, 
while for strongly localized solitons, with /i — > oo and s << 1, relation (j6|) demonstrates that Eq. (fT3|) yields 

e ( c ] ~ m/2, cf. Figs. land©. 

The general solution for the AS can be found after some algebraic manipulations with Eqs. © (replacing the 
equations by their sum and difference, and solving the latter equations for A 2 + B 2 and AB, cf. similar exact 
solutions obtained for the SSB in the double- well pseudopotential in Ref. (l2l|): 

A = £ y2.s(i - gj 
v / i + \/i — iPF' 

B = yi^^/l + y/i -4s 2 e 2 . (14) 

This solution exists at e < e c , where e c is precisely the bifurcation value given by Eq. (fF5]l . The conclusion that the 
asymmetric modes exist when the coupling constant (e) is not too large is very natural [13]- [11]) [HI- Indeed, the 
extreme case of the asymmetric solution, with A ^ and B = 0, is obviously possible in the limit of e = 0. Equally 
natural is the conclusion that the symmetric mode is stable at e > e c and unstable at e < s c , while the asymmetric 
mode is stable in the entire region of its existence, as seen in Fig. [2] Finally, we notice that the SSB bifurcations 
observed in Fig. [5] (and in Fig. [S] below, which pertains to System 2) are clearly the pitchfork bifurcations of the 
supercritical type, on the contrary to the subcritical bifurcation reported in Ref. [23 ] for the model based on two 
uniformly coupled parallel chains. 

For antisymmetric solitons (AnS) the amplitude is obtained from Eq. by substituting A = —B, which yields 

A 3 =(l-s*)(a- 1 +e). (15) 

Relation (|T3|) is plotted versus e for fixed /i = 5 by blue lines in Fig. [2] As follows from this relation, the VA predicts 
the existence of the AnS in the entire parameter space, on the contrary to the limited existence regions predicted 
for the SyS and AS modes. However, numerical results reported below demonstrate that the AnS also have their 
existence limits, see Fig. I^b). 



3. Stability 



The stability of the discrete solitons predicted by the VA can be first estimated by dint of the VK criterion, 
dPjd[i > 0, where P = X)n=-oo u n is the total power (norm) of the soliton [27| . Being necessary, but not the 
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sufficient, this condition should be combined with the spectral criterion, which requires all the eigenvalues of small 
perturbations around the solitons to be stable. 
The power corresponding to ansatz (j4|) is 

P=(l-s 2 y 1 (A 2 +B 2 ). (16) 

For the SyS solution with A — B, Eqs. (JIB} and (10) yield P = 2(1 + s 2 ) (s" 1 - e) . This expression satisfies condition 
dP/ds < 0, which is tantamount to the VK criterion, at all values of e. Nevertheless, it is obvious that the the 
branch of the symmetric solitons is destabilized by the SSB bifurcation, hence it is unstable at e < e f ,, as shown in 
Fig. [5] (it is well known that the respective instability is not detected by the VK criterion [13]- [HI; [Hi)- 

For the AS solutions, the substitution of Eq. (fT4")) into Eq. (fTB"|) produces a simple expression for the total power, 
which does not depend on e: 

P=(l + s 2 )/s. (17) 

It also satisfies the VK criterion in the entire region of the existence of the asymmetric interface modes. 

Finally, for the antisymmetric modes in System 1, with A = — B, the use of Eq. (|15l) gives P = 2(1 + s 2 ) (s _1 + e) . 
In this case, condition dP/ds < is satisfied only at 

s > £ « = (1 - ,s 2 ) S ~ 3 , (18) 

while the antisymmetric modes are VK-unstable in the other part of their existence region, at e < e s . This is in 
contrast with the above conclusions that the symmetric and asymmetric solitons comply with the VK criterion in 
their entire existence domains. 



4- Comparison to the uniform chain 

It is instructive to compare the results predicted for System 1 by the VA to the well-known properties of the uniform 
infinite chain described by the usual DNLSE, which corresponds to Eqs. (UJ with e = 1 [10]. To this end, one may 
consider cuts of Figs. H^a) and (b) through e = 1. Along these cuts, one finds an unstable SyS, stable AS, and stable 
AnS. In terms of the infinite uniform chain, the SyS corresponds to an inter- site- centered (alias off-site) discrete 
soliton, which is, indeed, always unstable in the usual DNLSE [l(| • The fact that the spring defect with e > 1 makes 
the off-site soliton stable beyond the bifurcation point, i.e., at e < e c , is remarkable by itself. 

Further, what we define as the AS, corresponds to the usual on-site-centered soliton in the ordinary DNLSE, which 
is always stable, as is well known. And finally, the AnS soliton is tantamount to the localized twisted mode in the 
DNLSE setting [2^, which is stable with small amplitudes and unstable if its amplitude is larger, as in the case 
corresponding to Fig. [^b) at s = 1. 

B. System 2: the parallel coupling 

The system formed by two identical infinite uniform chains, coupled by the transverse link at the single site, is 
shown in Fig. [4] 
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FIG. 4: The schematic of System 2, which is built of two uniform identical chains placed in parallel planes, which are linearly 
coupled in the transverse direction at site n = 0. The longitudinal and transverse coupling constants are C and e, respectively. 
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The corresponding DNLSE system is (where, as well as in Eqs. ([!}, we set C = 7 = 1, by means of rescaling) 



i-r^- + </Wi + 0n-i7 I0«| 2 4>n = 0, n ^ 0, 



+ 0i + 0-i +e-0o + |0o | 2 0o 



0, 



i-T^ +~4>n+l + V'n-l + = 0, n ^ 0, 

i— +^1 +'0-1 +£0o +7lV'o| 2 '0o = 0. 
az 



(19) 



The stationary solutions are looked for in the form of 0„ = w n exp(i/iz) and ip n = v n exp(ifj,z) , with a common 
propagation constant, //, and real functions u n and i> ra satisfying the following equations: 

- fiu n + u n+ i + u n -i + \u n \ 2 u n = 0; n ^ 

-fiv n + v n+1 +w„_i + \v n \ 2 v n = 0; n ^ 

— /iu + «i + u_i + ev + |"o| 2 "o = 0; 

-fivo + V\ + v-i + euq + \v Q \ 2 v q = 0. (20) 



1. The variational approximation 

The stationary equations can be derived from the Lagrangian, L = L\ + L2 + 2eu$vq, where L\ and Li are the 
intrinsic Lagrangians of the two uniform chains. Similar to System 1, the VA is based on the exponentially localized 
ansatz with possibly different amplitudes A and B but a common width, a" 1 , cf. Eqs. (j3]): 

u n = Acxp (— a\n\), v n — B exp (— a\n\). (21) 

By substituting this ansatz into the Lagrangian, and using, as above, Eq. ([6]) to eliminate fx, we obtain 

L cS = - — - (A 2 + 5 2 ) + ii±4 (A 4 + B 4 ) + 2eAB, (22) 

cf. Eqs. ([7J and (|SJ). Finally, the Euler-Lagrange equations, dL e g/dA = dL e s/dB = 0, applied to Lagrangian (j2"2"1) . 
take the following form: 

1 - s 2 1 + s 4 
s 1 — s 4 

— - — —B + ^——^-tB + eA = 0. (23) 
s 1 — s 4 



2. Existence regions for soliton solutions 
The solution for the SyS (symmetric soliton), with A = B, is immediately obtained from Eq. (|23|) : 

,2 1-S 4 /1-S 2 \ 

This relation is plotted in Fig. [5] for /.i = 5. Similar to the case of System 1, it follows from Eq. (|2~i)) that in System 2 

the SyS exists at s Kei 2 ^ (l — s 2 ) /s, as shown in Figs. [5]and|ni In the latter figure, the critical value is plotted 
by the black line. 

The existence range of the AS solutions can be determined by means of the same procedure which was used above 
in the case of System 1. The result is that, in the present case, the asymmetric modes exist at 

£ < 4 2) - (1 - * 2 ) / (2s) , (25) 
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cf. Eq. (p~3|) . For the broad asymmetric solitons (s — > 1), approaches the limiting value e^= 0, while for the 
strongly localized solitons (s << 1) we obtain e c ps /x/2. The full solution for the AS mode can be found too, cf. Eq. 
([13]) for System 1: 



A = 



B = e 



1 



2s(l + s 4 ) 
2s(l - s 4 ) 



(1 - s 2 ) + V(l - s 2 ) 2 - 4s 



2 C 2 



(26) 



4s 2 e 2 



Note that the AS solution exists at e < (l — s 2 ) / (2s), which is exactly equivalent to Eq. (|25j) . obtained above by 
means of a different algebra. The amplitudes of AS at the interface sites are plotted in Fig. [5] 
The solution for the antisymmetric (A = —B) solitons (AnS) produced by Eqs. (f23| is 



1 - s 4 / 1 - s 



(27) 



As seen from here, the AnS are predicted by the VA to exist in the entire parameter space. 
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FIG. 5: (Color online) The same as in Fig. [2jb) (i.e., with /i = 5), but for System 2. 
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3. Stability 

The VK criterion may also be applied to the soliton solutions predicted in System 2. The total power of the 
solutions based on ansatz (|2"TT) is 

P=\±^{A 2 +B 2 ). (28) 
1 — s l 

For the SyS solution ([MD with A = B, Eq. {25]) yields 



1 + s 4 

which satisfies the VK condition, dP/ds < (tantamount to dP/dfi > 0) in the entire existence region. For the AS 
solution given by Eq. (|26p . the calculation of norm (f2"5|) yields a result which does not depend on e, cf. Eq. p7|) for 
the AS solution in System 1: 

P ^ 1 + s2 ^-^ (30) 
F ~ s (l + ,s 4 ) ' 1 j 

whose VK slope, dP/d^i, is again positive in the entire region where this solution exists. 
Lastly, for the AnS solution, with A — —B, the power (|29[) is 



p = H(l±i!)!(i^! +l| , (31) 



1 + s 4 

whose slope may change its sign. The eventual conclusion is that the AnS solution satisfies the VK criterion only in 
a part of its existence region, viz., at 

£ > £ (2) _ 1 + g2 _ izi! (32) 

£>£s "4 S 3(1_ S 2 } , > M 

cf. similar result (fT8|) obtained above for the AnS solutions in System 1. 



III. NUMERICAL RESULTS 



To verify the predictions of the VA, we numerically solved stationary equations ([2]) and (|20|) for both models, 
applying the relaxation method and an algorithm based on the modified Powell minimization method [29]. The 
initial guess used to construct fundamental solitons centered at the interface of the linked semi-infinite chains which 
constitute System 1 (Fig. [Tj) was taken as U Q = U\ — A > for SyS solutions, U = A > 0, U\ = B > for AS 
solutions, and Uq = A > 0, U% — — A for ones of the AnS type, with the VA-predicted values of A and B, while in 
initial values of the lattice field are set to be zero at all other sites. 

The initial ansatz for the parallel-coupled lattices which constitute System 2 (Fig. [4j was Uq = Vq = A > 
for SyS modes, U a = A > 0, V Q = B > for AS modes, U = A > 0, V Q = -A for modes of the AnS type, and 
U n = V n = at all other sites. The results presented here are obtained for identical coupled chains with Ni = N2 = 50 
or Ni = N2 = 51 sites, for Systems 1 and 2, respectively. In other words, the total number of sites is N — 100 in 
System 1, and N = 102 in System 2, which features the parallel chains. The link between the chains in System 1 was 
set between the sites with indexes n = and n = 1, and in System 2 the transverse coupling between the parallel 
lattices was introduced at n = 0. 

The stability of the stationary modes was first checked by dint of the linear-stability analysis, i.e., the calculation of 
eigenvalues (EVs) for modes of small perturbations, following the lines of Refs. [2!|. The respective calculations were 
performed in parameter space (e, fj,). Then, the results were verified in direct numerical simulations of full equations, 
([T]) and (I19[) . The simulations were based on a numerical code which used the sixth-order Runge-Kutta algorithm, as 
in Ref. [29(. The simulations were initialized by taking the stationary soliton profiles, to which random perturbations 
were added. 

Typical shapes of symmetric, asymmetric and antisymmetric solitons found in the numerical form are displayed in 
Fig. [7] for System 1 and in Fig. [S]for System 2. In the same figures, the numerical shapes are compared with those 
obtained by means of the VA. The respective dependencies of the solitons' amplitudes A and B on parameter e, which 
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FIG. 7: Profiles of the fundamental solitons (interface modes) in System 1: (a) symmetric (SyS), (b) asymmetric (AS), and (c) 
antisymmetric (AnS). The inter-chain linkage strength is e = 0.65. The total powers corresponding to each mode are indicated 
in the panels. Black solid lines with squares denote numerically generated solitons, while red dashed lines with circles denote 
variationally obtained solitons. 
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FIG. 8: (Color online) Profiles of fundamental solitons in System 2 (two parallel lattices, Li and L2, coupled by the transverse 
link with strength e), obtained for /1 = 6. The values of e and the total powers are indicated in the panels. Black solid and 
red dashed lines show the profiles in lattices Li and L2 obtained numerically. The corresponding soliton profiles predicted by 
the VA are denoted by hollow squares (green for Li and blue for L2). (a) A symmetric (SyS) soliton, (b) an asymmetric (AS) 
soliton, and (c) an antisymmetric (AnS) soliton. 



accounts for the linkage between the two chains, are displayed for all the soliton species, alongside their VA-predicted 
counterparts, in the above figures [2] and [3J Further, the numerical results demonstrate that all the species of the 
solitons, including the AnS, exist in bounded regions of the parameter space, as shown in the above figures [3] and |H1 
for Systems 1 and 2, respectively. 

The comparison with the numerical results demonstrates that the predictions of the VA for the symmetric and 
asymmetric discrete solitons are very accurate. However, the VA fails to predict the existence borders for the anti- 
symmetric modes. A reason for the latter problem may be that the strong interaction force at the interface, in the 
case of the opposite signs of the lattice fields at adjacent sites, see Figs. [7£c) andlHJc), makes the simple form of the 
VA adopted above inaccurate. 

In the instability region, the SyS solutions are characterized by the EV spectrum consisting of real pairs. Under 
small asymmetric perturbations, an unstable SyS sheds off a part of its total power and relaxes into an asymmetric 
breather with a smaller power, as shown for System 1 in Fig. [9] (a). As in other conservative systems, 
the asymmetric breathers cannot readily self-trap into stationary solitons. In System 2, the behavior of the unstable 
symmetric solitons is quite similar. 

In both Systems 1 and 2, the EV spectra for the stationary asymmetric modes are stable in the entire existence 
region. This is also confirmed by direct simulations, which show that perturbed AS develop small-amplitude persistent 
oscillations without any trend to destruction. 

The antisymmetric solitons change their stability twice, in both configurations considered, 1 and 2. With the 
decrease of e at fixed /u, the unstable AnS branch, which is characterized by a pair of real EVs [in the area between 
the green dotted and black solid lines in Figs. [3j(b) and|6](b)], changes into the stable one. With the further decrease 
of e, the AnS loses its stability through the appearance of two pairs (a quartet) of complex EVs with significant 
real parts [in the area below the red dashed line in Figs. [3] (b) and|5] (b)]. For example, the antisymmetric modes 
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FIG. 9: Typical examples of the evolution of perturbed unstable solitons in System 1. (a) A symmetric soliton, with e = 
0.65, /i = 6.425. Other panels display antisymmetric solitons, with fi = 6: (b) e = 9.1 (real unstable EVs, see Fig. [5J, and (c) 
e = 1.15 (complex EVs). The unstable symmetric soliton, after shedding off a part of its power, evolves into an asymmetric 
breather. The unstable antisymmetric solitons transform into breathers too, see the text. 



corresponding to fx = 5 acquire the oscillatory instability, with the decrease of e, at e 1.5 and 2 in Systems 1 and 
2, respectively, which is indicated by the appearance of a quartet of complex EVs in their spectra (in the areas below 
the red line in Figs [3] and [5]) . Note that the VA has also shown that the AnS change their stability according to the 
VK criterion, see Eqs. (TT51) and (|32l) . However, the corresponding VA-predicted critical curves, ei 1 ' 2 ^^), are found 
to be situated beyond the numerically obtained existence regions of the AnS solutions. 

The predictions concerning the stability of the AnS were checked by direct simulations. In Fig. IHJb) and (c), 
the evolution of typical unstable antisymmetric modes in System 1 is displayed. Panels (b) and (c) in this figure 
illustrate the evolution of the unstable AnS whose EV spectrum contains a pair of real EVs, or a complex EV quartet, 
respectively. The unstable antisymmetric modes radiate away a part of their power, relaxing to antisymmetric (b) or 
asymmetric (c) interface modes, which exhibit small-amplitude persistent oscillations. In System 2, the same happens 
with the unstable AnS belonging to the instability area below the red (dashed) line in Fig. [6] However, in contrast to 
that, the AnS in System 2 which is characterized by the EV spectrum with a pair of real EVs turn out to be robust 
under the action of small perturbations. They do not emit radiation waves, staying localized and strongly pinned to 
the link connecting the two infinite chains. 

Returning to the global existence diagrams, it is worthy to note that two bistability areas can be identified in both 
systems, 1 and 2: the domain of the coexistence of stable symmetric and antisymmetric solitons, or the one featuring 
the simultaneous stability of asymmetric and antisymmetric modes, on the opposite sides on the SSB bifurcation. 
This result is in accordance with similar findings reported in other linearly-coupled two-component systems featuring 
the self- focusing nonlinearity [111, [22[ • 



IV. CONCLUSION 



In this paper, we have investigated several species of fundamental localized modes formed at the interface between 
two linearly coupled lattices (chains) with the cubic on-site nonlinearity. Two configurations were considered: the 
one with the linkage between two semi-infinite chains (in other words, the usual discrete NLS model with a spring 
defect), and two infinite chains placed in parallel planes which are coupled by the transverse link at one site. These 
systems can be readily implemented as arrays of optical waveguides. 

In both models, the VA (variational approximation) was used to predict the existence regions, in the parameter space 
of (n,e) (the propagation constant and the strength of the coupling between the two subsystems), for the localized 
symmetric, asymmetric and antisymmetric discrete solitons pinned to the interface. The stability of the modes was 
predicted as per the VK criterion and general properties of the supercritical pitchfork bifurcation, which destabilizes 
the SyS (symmetric solitons) giving birth to AS (asymmetric solitons) . The predictions were verified against numerical 
results, as concerns the existence and stability of all the soliton species. In both systems considered, the existence 
regions of all the localized modes are bounded. The AS are stable in their entire existence domain. The existence 
and stability domains for AnS (antisymmetric solitons) were found in the numerical form. Both systems give rise to 
the bistability between the AnS on the one hand, and either SyS or AS, on the other. Direct simulations demonstrate 
that unstable SyS are transformed into asymmetric breathers. Those antisymmetric modes which are unstable shed 
off a part of the total power and also evolve into breathers. In System 2, in a part of the parametric domain where 
the computation of the eigenvalues predicts the instability of the antisymmetric localized modes, they actually evolve 
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into strongly pinned robust spikes. 
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